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Abstract 

We  introduce  a  frequentistic  validation  framework  for  assessment  —  acceptance  or  rejection  —  of  the  con¬ 
sistency  of  a  proposed  parametrized  partial  differential  equation  model  with  respect  to  (noisy)  experimental 
data  from  a  physical  system.  Our  method  builds  upon  the  Hotelling  T2  statistical  hypothesis  test  for  bias 
first  introduced  by  Bald  &  Sargent  in  1984  and  subsequently  extended  by  McFarland  &  Mahadevan  2008. 
Our  approach  introduces  two  new  elements:  a  spectral  representation  of  the  misfit  which  reduces  the  di¬ 
mensionality  and  variance  of  the  underlying  multivariate  Gaussian  but  without  introduction  of  the  usual 
regression  assumptions;  a  certified  (verified)  reduced  basis  approximation  —  reduced  order  model  —  which 
greatly  accelerates  computational  performance  but  without  any  loss  of  rigor  relative  to  the  full  (finite  el¬ 
ement)  discretization.  We  illustrate  our  approach  with  examples  from  heat  transfer  and  acoustics,  both 
based  on  synthetic  data.  We  demonstrate  that  we  can  efficiently  identify  possibility  regions  that  characterize 
parameter  uncertainty;  furthermore,  in  the  case  that  the  possibility  region  is  empty,  we  can  deduce  the 
presence  of  “unmodeled  physics”  such  as  cracks  or  heterogeneities. 

Keywords:  Validation;  uncertainty  quantification;  hypothesis  testing;  partial  differential  equation;  reduced 
basis  method;  a  posteriori  error  bound 


1.  Introduction 

Parameterized  models  play  a  central  role  in  engineering.  Our  focus  here  is  parametrized  partial  differ¬ 
ential  equations  (PDEs)  that  arise  for  example  in  heat  transfer,  solid  mechanics,  acoustics,  fluid  flow,  and 
electromagnetics.  The  parameters  appear  as  coefficients  in  the  partial  differential  equation  and  may  represent 
physical  properties,  interactions  with  the  environment,  or  geometry.  Often  one  set  of  parameters  will  serve 
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to  relate  the  PDE  model  to  the  actual  Physical  System  (PS)  of  interest,  while  another  set  of  parameters  will 
then  serve  to  subsequently  design  or  optimize  or  control  the  Physical  System.  Our  interest  in  this  paper  is 
the  former. 

There  are  many  frameworks  by  which  to  relate  a  proposed  PDE  model  to  the  Physical  System  of  interest. 
Parameter  estimation  [1,  2]  and  calibration  (or  “update”)  [3,  4,  5]  search  for  a  best  choice  of  parameters 
relative  to  (typically  noisy)  experimental  data.  Validation  [3]  assesses  the  consistency  of  any  candidate 
parameter  value,  and  hence  associated  proposed  PDE  model,  with  (typically  noisy)  experimental  data.  Both 
activities  —  calibration  and  in  particular  validation  —  play  a  central  role  in  Uncertainty  Quantification 
(“UQ”).  Our  focus  in  this  paper  is  on  validation. 

Validation  procedures  fall  into  two  broad  categories:  Bayesian  and  frequentistic  [6].  Bayesian  validation 
builds  directly  upon  Bayesian  parameter  estimation  [2,  4]  —  from  experimental  data  to  likelihood  through 
prior  to  posterior  and  finally  credible  regions  [6].  The  Bayesian  approach  is  quite  general  and  can  yield 
actionable  information  on  many  quantities  of  interest;  however,  the  prior  can  be  difficult  to  specify  or 
interpret  [3,  7],  and  hence  absolute  conclusions  can  remain  elusive.  Frequentistic  validation  procedures  [8,  7] 
are  typically  constructed  as  statistical  tests  of  appropriate  hypotheses  on  bias  which  then  permit  acceptance 
or  rejection  of  a  particular  choice  of  candidate  parameters.  The  frequentistic  approach  avoids  the  Bayesian 
prior  and  can  thus  provide  more  rigorous  or  at  least  less  subjective  conclusions;  however,  Type  II  errors  [9] 
can  typically  not  be  controlled,  and  hence  the  implications  are  perforce  restricted. 

In  this  paper  we  consider  frequentistic  validation.  The  seminal  work  in  frequentistic  validation,  due 
to  Bald  and  Sargent  [8],  is  based  on  the  multivariate  Gaussian  Hotelling  T2  test  [10].  More  recently, 
McFarland  and  Mahadevan  [7]  extend  and  apply  the  Bald  &  Sargent  approach  to  the  thermal  validation 
challenge  problem  proposed  by  Sandia  National  Laboratory  [11].  We  shall  consider  in  this  paper  two  further 
innovations  to  the  methods  proposed  in  [8,  7]:  the  first,  (arguably)  an  improvement,  is  introduced  to  reduce 
the  dimensionality  of  multivariate  Gaussian;  the  second,  an  extension,  is  introduced  to  permit  more  efficient 
computation.  We  discuss  each  in  turn. 

As  regards  the  first  innovation,  we  develop  the  hypothesis  test  based  not  on  pointwise  data  but  rather  on 
“spectral”  Legendre  coefficients  [12].  This  approach  has  two  advantages:  it  leads  to  a  significant  reduction 
in  the  dimensionality  of  the  multivariate  Gaussian  and  hence  permits  more  rigorous  inference  in  particular 
of  the  variance  (e.g.,  without  further  hypotheses);  it  permits  inferences  from  fewer  experimental  realizations 
(or  replications)  since  the  spectral  coefficients  implicitly  average  the  data  and  hence  reduce  the  variance,  at 
least  for  correlation  scales  small  compared  to  the  full  observation  window. 
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We  emphasize  that  our  approach  differs  from  standard  regression  [13]  in  that  we  do  not  assume  that 
our  truncated  spectral  expansion  is  unbiased:  our  inferences  derive  from  a  Linear  Optimization  problem 
constrained  by  the  Legendre  coefficient  joint  confidence  intervals.  Note  also  that  we  do  not  assume  uncor¬ 
related  noise  and  in  fact  we  permit  any  Gaussian  process  [14]:  of  course,  unlike  the  standard  regression 
framework,  we  can  not  simultaneously  deduce  bias  and  variance  from  a  single  realization;  nevertheless,  as 
already  indicated,  relatively  few  replications  suffice. 

Validation  perforce  requires  the  consideration  of  many  candidate  parameter  values.  In  the  PDE  context, 
the  evaluation  of  the  model  by  classical  (in  our  case)  finite  element  [15]  methods  can  be  very  expensive.  Our 
second  innovation  is  thus  to  replace  the  “truth”  finite  element  discretization  with  a  less  costly  surrogate:  a 
certified  reduced  basis  approximation  [16,  17,  18,  19,  20,  21]  which  provides  not  just  an  inexpensive  output 
prediction  but  also  —  in  the  sense  of  “verification”  -  an  inexpensive  but  rigorous  output  error  bound. 
The  reduced  basis  (RB)  approach  is  efficient  in  the  many-query  and  real-time  contexts,  both  of  which 
are  relevant  in  (calibration  and)  validation  procedures.  In  order  to  incorporate  the  certified  reduced  basis 
method  into  the  frequentistic  validation  framework  we  consider  a  composite  hypothesis  [9,  8].  The  RB  error 
bound  is  crucial  in  ensuring  proper  distinction  between  experimental  noise,  bias  in  the  model  prediction, 
and  numerical  (reduced  basis)  error:  without  the  RB  error  bound,  we  might  conflate  numerical  error  with 
bias  and  hence  reject  parameter  values  which  in  fact  are  consistent  with  experimental  data.  (We  note  that 
metamodels  and  reduced  models  are  gainfully  invoked  [22,  7,  5,  23]  in  several  calibration  and  validation 
approaches,  but  only  in  a  few  studies  [24]  are  rigorous  error  bounds  considered.) 

The  rest  of  this  paper  is  arranged  as  follows.  In  Section  2  we  introduce  the  statistical  framework 
a  hypothesis  test  for  bias  -  for  our  approach  independent  of  any  specific  mathematical  model  or  physical 
system.  In  Section  3  we  specialize  the  framework  to  the  case  in  which  our  model  is  given  by  a  proposed 
parametrized  PDE;  in  particular,  we  incorporate  the  reduced  basis  output  and  reduced  basis  output  error 
bound  into  our  hypothesis  test.  In  Section  4,  we  present  numerical  results  for  a  transient  thermal  conduction 
problem  and  a  Helmholtz  acoustics  problem.  In  both  cases  we  demonstrate  the  ability  of  our  approach  to 
identify  “possibility  regions”  based  on  many-query  parameter  sweeps  over  the  parameter  domain,  and  also 
to  deduce  the  presence  of  “unmodeled  physics”  (defects)  in  a  Physical  System. 

2.  A  Statistical  Result 

We  introduce  a  compact  domain  Dv  £  ;  in  this  paper  we  shall  exclusively  consider  P„  =  1,  and  hence 

we  may  write  T>v  =  [r'min ,  r'max]  •  We  then  define  a  (centered)  Gaussian  random  process  [14]  G(v,u>)  such 
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that  the  mean  is  zero, 


J  dP(u)G(v,uj)  =  0,Vv  e 


and  the  average  variance  is  unity, 


^max  ^min 


—  f  c?P(w)  f  G2(y,  u>)dv  =  1. 
in  J  JVU 


Finally,  for  given  positive  definite  covariance  kernel  C,  G  satisfies 


J  dP(w)G(v,  w)G(j/,w)  =  C(v,  z/), 


(1) 


(2) 


(3) 


Here  tflP(w)  indicates  integration  with  respect  to  the  probability  measure.  Note  that  the  variance  condition 
on  G,  (2),  represents  a  scaling  constraint  on  the  covariance  C. 

In  actual  practice  such  a  Gaussian  process  can  be  represented  by  a  Karhunen-Loeve  expansion  [25].  We 
first  find  \  '■  TV  — t  and  A  £  K  solutions  of  the  eigenproblem 


C{v,  v')x(v')dv' 


Ax(^) 


(4) 


with  normalization  fv  X1  {v)dv / {vmax  —  ;n)  =  1.  We  may  then  express  our  random  process  as 


G(u,u)  =  '/^iWi{u)xi(v),  (5) 

i=i 

where  the  Wi  are  independent  normally  distributed  random  variables  with  zero  mean  and  unity  variance. 
Note  we  order  the  (positive)  eigenvalues  in  decreasing  magnitude  (Ai  >  A2  •  •  •  )•  The  sum  of  the  eigenvalues 
is  the  integrated  variance,  and  thus  from  (2)  we  conclude  that 


1 

^max  ^min 


Xi  =  L 
1=1 


In  practice  we  truncate  the  expansion  (5)  at  L  terms  such  that  sum  of  the  remaining  eigenvalues  Al+i,  \l+2,  ■  ■  ■ 
is  small  compared  to  unity. 

In  our  procedure  we  shall  assume  that  the  underlying  noise  or  stochastic  element  is  a  Gaussian  process. 
However  —  and  unlike  kriging  [26]  —  we  will  not  need  to  know  the  details  of  the  process  in  order  to 

construct  our  approximation  or  provide  our  statistical  inferences.  Nevertheless,  for  purposes  of  generating 


4 


representative  synthetic  data,  we  will  consider  two  convenient  correlation  functions: 

(i)  White  noise: 

{Const,  v  =  v' 

(6) 

0,  v^v' 

(ii)  Correlated  noise: 

C(is ,  is')  =  Const  exp(— 1^  —  is'  |/A„),  (7) 

where  Av  is  a  specified  correlation  scale  in  v. 

In  each  case  Const  is  chosen  to  ensure  unity  average  variance  over  T>v,  (2).  With  increasing  correlation  scale 
the  spectrum  of  (4)  decreases  more  rapidly  and  we  may  choose  L  correspondingly  smaller  —  we  may  include 
fewer  and  fewer  terms  in  the  truncated  Karhunen-Loeve  expansion. 

We  now  introduce  a  function  S  :  T>v  — >  R.  We  then  define,  for  “given”  integrated  variance  e,  the  random 
process 

Y(is,  w)  =  6  (is)  +  eG(is,  w).  (8) 

Note  that  the  expectation  of  Y(is,  •)  is  6 (is)  from  our  assumption  (1)  on  G.  We  emphasize  that  in  our  inference 
procedure  we  will  not  need  to  know  the  magnitude  of  the  average  (over  D„)  variance,  e. 

We  next  introduce  a  set  of  “quadrature  points”  (ultimately  measurement  points)  isrn  g  Vv.  1  <  m  <  M, 
and  associated  positive  quadrature  weights  pm,  1  <  in  <  M.  We  then  provide  a  “projection  matrix” 

w®  = - - U  (-1  +  2  (  \\  .  1  <  m  <  M,  0  <i  <  I,  (9) 

^max  ^min  \  \  ^max  ^min  /  / 

where  Li  :  [—1,1]  — »•  R  is  the  Legendre  polynomial  of  order  i  [27].  The  “Legendre  coefficients”  of  6  are  then 
given  by 

M 

=  ^2  WmS(Vrn),  0  <  *  <  I.  (10) 

m=  1 

We  note  that  if  6  is  a  smooth  function  then  we  expect  rapid  decrease  in  the  Legendre  coefficients. 

There  is  much  freedom  in  the  choice  of  quadrature  (points  and  hence  weights).  In  this  paper  we  employ 
the  trapezoidal  rule: 

=  ^min  +  (m  -  1)A IS,  pm  =  CmAlS ,  1  <  TO  <  M,  (11) 

where  A  is  =  (vmax  —  ismin)/(M  —  1),  and  c\  =  Cm  =  1/2,  cm  =  1  for  2  <  m  <  M  —  1.  Clearly  there  may  be 
some  advantage  to  a  Gauss-Legendre  (or  Gauss-Lobatto)  quadrature  [27],  however  the  associated  clustered 
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point  distribution  could  be  unnatural  in  an  experimental  implementation  and  we  thus  restrict  attention  to 
low-order  schemes. 

Now  let 

Y^P(co)  =  Y(um,co),  1  <m<M,  (12) 


be  a  realization  ui  of  our  random  process  (8):  a  particular  set  of  experimental  measurements  of  Y(i/,u>)  at 
the  points  vmi  1  <  m  <  M.  We  may  then  introduce 

M 

Z[iIM  =  £  w^Y^P M,  0  <  t  <  /;  (13) 

m—1 

we  may  view  as  a  (noisy)  approximation  to  the  Legendre  coefficient  <5^  of  (10).  It  is  instructive  to 

rewrite  (13)  in  terms  of  our  Karhunen-Loeve  expansion  (now  truncated  at  L  terms):  from  (5),  (8),  and  (13) 
we  obtain 

L  M 

zW(u)=6M+eJ2Wi(u)[^Y,™[£xi(vm)],  0  <i<I.  (14) 

l—l  m= 1 

It  then  follows  that  the  Legendre  coefficient  random  vector,  {Z^°\  Z I1! , . . . ,  Z^)(ui),  is  (/+l)-variate  normally 
distributed  with  mean  (d ^ , . . . ,  )  [9] . 

We  next  introduce  an  ensemble  of  I\  realizations  (presumed  of  course  independent)  oj\, ,  u>k  '■  in  practice 
this  ensemble  will  take  the  form  of  K  sets  of  M  measurements  each.  We  may  then  introduce  the  sample 
mean  of  the  Legendre  coeffients  over  this  ensemble  as 

1  K 

Zk  =  xJ2Z[i]^'  0<i<I,  (15) 

k=l 

for  Z^(u>)  given  by  (13).  Similarly,  the  sample  variance  over  the  ensemble  is  given  by 

=  TE(Z[l'W-^)2-  0  <i<I.  (16) 

k— 1 


(We  shall  not  need  the  full  sample  covariance  for  reasons  described  below.) 
It  then  follows  that  the  intervals 


T 


K 


0  <i<I, 


(17) 
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represents  a  rectangular  confidence  region  which  encloses  the  strict  elliptical  confidence  region  [10]  —  the 
ith  dimension  of  this  rectangular  region  is  derived  by  maximizing  the  ith  component  of  y  £  ]R7+1,  for 
i  =  1, ...,/  +  1,  with  y  constrained  to  reside  on  the  elliptical  confidence  region  from  [10].  We  prefer  the 
rectangular  region  in  anticipation  of  the  Linear  Program  introduced  below;  ultimately,  however,  the  elliptical 
region  (or  an  oriented  rectangular  region)  could  be  incorporated  and  would  of  course  sharpen  our  results. 
Second,  we  emphasize  that  our  inference  related  to  the  I  +  1  Legendre  coefficients  does  not  assume  -  in 
contrast  to  a  corresponding  (I  +  l)-term  regression  model  [13]  -  that  the  bias  is  zero:  our  inference  (18) 

does  not  suppose  that  5(is)  is  exactly  represented  by  the  first  ( I  +  1)  terms  of  a  Legendre  expansion.  Third, 
we  observe  that  typically  the  number  of  measurements  M  might  be  quite  large  compared  to  the  requisite 
number  of  Legendre  coefficients  I  +  1;  this  will,  in  view  of  the  denominator  in  (19),  permit  inference  with 
fewer  replications  than  in  pointwise  approaches  (e.g.,  McFarland  &  Mahadevan  [7]).  Finally,  we  implicitly 
obtain  variance  reduction  since  the  Z M  are  in  some  sense  “moments”  constructed  as  averages  over  the  M 
measurements;  this  further  reduces  the  number  of  replications  required. 

We  elaborate  further  on  the  variance  reduction.  We  recall  that  the  procedure  defined  above  requires  an 
ensemble  of  K  realizations  (which  we  may  also  view  as  repetitions  or  replications)  each  of  which  comprises  M 
measurements.  In  general,  given  that  we  allow  correlated  errors  —  often  the  case,  for  example  if  v  represents 
time  or  space  and  disturbances  and  measurement  errors  are  spectrally  broad  —  a  single  realization  can  not 
suffice.  As  a  simple  counterexample  we  might  consider  the  case  in  which  G  has  infinite  correlation  scale:  a 
single  realization  can  not  possibly  distinguish  signal  from  noise.  However,  our  procedure  does  nevertheless 
benefit  from  the  M  data  points  associated  with  each  realization:  if  G  has  relatively  short  correlation  scale 


7 


A„  —  in  the  limit  A„  — >  0  we  recover  the  standard  regression  assumption  [13]  of  white  Gaussian  noise  - 
then  the  variance  of  the  pseudo  sample  means  (w),  0  <  i  <  I,  will  be  small  and  only  very  few  realizations 
will  be  required  in  order  to  obtain  an  acceptably  tight  confidence  interval. 

We  now  define 

^max  =  max  |(^(z/m)|-  (20) 

We  shall  then  wish  to  test  the  hypothesis 

n  =  { <  b  }  (21) 

for  some  positive  constant  B.  We  may  view  (21)  as  a  hypothesis  on  bias;  we  introduce  in  the  next  section 
the  particular  interpretation  of  bias  relevant  in  our  context. 

We  now  describe  the  rejection  criterion.  We  shall  reject  hypothesis  TL  if  and  only  if 

Lx  >  B ,  (22) 


where 

^max  =  min  max  \vm\,  (23) 


and  the  constraint  set  7Z  is  given  by 


M 

n  =  {v  G  rm  I  Y,  w™v™  6ll'0<j<  I}-  (24) 

m—1 


We  shall  denote  by  1  —  7  (where  7  here  determines  T7  in  (19))  the  size  of  the  test. 

In  particular,  we  now  derive  a  bound  for  the  Type  I  [9]  error  in 

Proposition  2.1.  The  probability  of  rejection  of  hypothesis  TL  when  hypothesis  TL  is  in  fact  true  is  less  than 
or  equal  to  1  —  7. 

Proof.  We  consider,  as  dictated  by  our  hypothesis,  the  case  in  which  the  hypothesis  TL  is  rejected:  jmax  >  B. 
We  first  suppose  that  the  condition  £  =  {<5 1*1  €  <  i  <  1}  is  satisfied,  which  irnples  that  S(um),  1  < 

m  <  M,  is  in  the  admissible  set  1Z.  Since  ^max  is  defined  in  (23),  (24),  as  the  minimum  over  v  €  1Z  of 
rnaxme{1;  ,jM}  \vm\,  it  follows  that  <5max  >  5max. 

Then,  under  the  above  conditions,  we  may  conclude  that  i5max  >  B ,  which  implies  that  the  hypothesis 


TL  is  indeed  false.  Hence  £  true  implies  correct  rejection  of  H,  and  thus  a  necessary  condition  for  incorrect 
rejection  of  hypothesis  TL  (Type  I  error)  is  £  false.  We  conclude  from  (18)  that  the  probability  of  a  Type  I 
error  is  less  than  or  equal  to  1  —  7.  □ 

Remark  2.1.  We  can  not  control  the  Type  II  error  [9]  —  in  which  we  accept  hypothesis  TL  when  hypothesis 
TL  is  in  fact  false  —  though  we  do  anticipate  that  for  larger  K  (tighter  confidence  intervals)  and  larger  I 
(better  estimates  for  Smax)  the  probability  of  Type  II  errors  shall  decrease.  Similarly,  we  are  not  able  to 
quantify  the  power  of  our  test. 

Finally,  we  note  that  our  test  (22)— (24)  can  be  efficiently  computed.  It  is  a  standard  procedure  to 
reformulate  (23),  (24),  as  a  simple  Linear  Program  (LP)  in  M  +  1  variables  and  2 (M  +  I  +  1)  constraints  as 
follows:  <5max  =  mine  over  (c,v  1,  •  •  • ,  Vm)  G  Rm+1,  where  (c,v  1, . . . ,  vm)  G  Rm+1  satisfies 

—c<vm<c,  m  =  1, ... ,  M, 

T[i]  /  V'  Id  /  -rid  •  _  n  T 

^K, min  —  WmVm  <  I^maxj  l  —  0,  ...  ,1. 

m—  1 

Note  each  two-sided  inequality  represents  two  constraints. 

3.  PDE-Based  Statistical  Inference 

We  now  proceed  to  the  parametrized  partial  differential  equation  context.  In  addition  to  v  £  T>v  which 
is  a  control  parameter  (or  set  of  parameters)  we  also  introduce  a  parameter  k  £  T>K  C  Rp"  (ultimately,  we 
shall  wish  to  determine  n  from  data).  We  are  given  a  proposed  parametrized  partial  differential  equation 
model  A4!',K;  this  partial  differential  equation  yields  an  exact  (model- prediction)  output  s:  Vu  x  T>K  — >•  R. 
In  actual  practice  s  may  not  be  obtained  analytically,  and  instead  we  resort  to  a  high-fidelity  finite  element 
[15]  “truth”  approximation  (here  h  denotes  the  mesh  size);  this  approximation  yields  a  truth  output 

Sh  '■  x  VK  — >  R.  We  presume  that  the  truth  approximation  is  sufficiently  rich  —  sufficiently  many  and 
judiciously  chosen  degrees  of  freedom  (. h  sufficiently  small)  —  such  that  the  difference  between  s  and  sh  is 
negligible.  Unfortunately,  this  truth  approximation  will  also  often  be  prohibitively  costly:  each  evaluation 
(u,  k)  — >  Sh(v,  k)  will  be  very  expensive. 

We  thus  introduce  a  certified  reduced  basis  approximation  [20,  28,  29,  30]  which  takes  advantage  of  the 
parametric  definition  and  parametric  structure  of  the  problem  to  greatly  reduce  the  computational  cost  in  the 
many-query  and  real-time  contexts.  We  shall  denote  the  N  degree-of-freedom  reduced  basis  approximation 
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by  M^K  ;  this  reduced  basis  approximation  yields  a  reduced  basis  output  sjv  :  XV  x  XV  — >  M.  and  a  reduced 
basis  output  error  bound  An  :  XV  x  XV  — >  ffi.  such  that  \sh{v,  At)  —  sn{v,  ac)|  <  Ajv(g  k),  V(i/,  At)  G  XV  x  XV- 
The  reduced  basis  method  accommodates  an  Offline-Online  computational  strategy  which  accepts  “many- 
truth”  cost  in  the  Offline  stage  in  exchange  for  very  rapid  evaluation  (v,k)  —a  Sn{v,  At),  An(v,  At)  in  the 
Online  (many-query,  real-time)  stage.  Note  that  in  actual  practice  we  consider  a  hierarchy  of  reduced  basis 
approximations  ,  1  <  IV  <  Nmax. 

We  assume  that  our  experimental  output  is  given  by  the  random  Gaussian  process 

V(v,  oj)  =  sPS(i/)  +  eG(u,  u>)  (25) 

where  sPS(^)  represents  the  response  of  our  physical  system.  Note  that  e  and  G  are  unknown  to  us:  we 
presume  only  that  the  random  process  G  is  Gaussian. 

We  now  wish  to  compare  the  output  of  our  proposed  parametrized  PDE  model  to  the  measured 

data  V(v,uj).  We  thus  form  the  new  random  variable  D(y,n,u>)  =  V{v,ui)  —  Sn(v,k)  or 

D{y7  k,lo)  =  sPS(iy)  —  Sn(v7  k)  +  eG(v,uj).  (26) 

We  then  presume  that  we  are  given  some  candidate  value  ac  =  k  in  VK  and  define  YK(y,u>)  =  D{v1  k)  or 

YK{v,  oj)  =  [sPS(ia)  -  Sn{v,  k)]  +  eG{v,u).  (27) 

We  note  that  (27)  has  the  same  form  as  (8)  but  now  S(f)  =  SK( v)  =  sPS(i/)  —  sn{v,k)  has  the  particular 
interpretation  of  the  deviation  or  misfit  between  experiment  and  (reduced  basis)  prediction. 

We  now  proceed  to  choose  the  number  of  measurements  per  realization,  M,  the  associated  measurement 
and  quadrature  points  um,pm,l  <  m  <  M,  the  order  of  the  Legendre  approximation  I,  and  finally  the 
sample  size  K.  We  then  perform  the  necessary  experiments  in  order  to  evaluate  the  zW(wfc),  1  <  k  < 
K,  0  <  i  <  I,  of  (13),  and  subsequently  the  Z^.  0  <  i  <  I.  of  (15).  Note  that  in  this  context  (12) 
should  be  interpreted  (for  given  k)  as  Y^xp(w)  =  Vexp(i zm,w)  —  SN{ym,k):  for  1  <  m  <  M ,  we  measure 
Vexp(ism,oj)  =  sPS(vm)  +  eG(vm,u))  and  then  subtract  the  reduced  basis  prediction  SN{ym,k). 

Finally,  we  now  set 

BK  =  max  Aj (28) 
and  we  consider  the  hypothesis  (21)  with  B  =  BK:  the  bound  in  our  test  is  now  related  to  the  maximum 
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reduced  basis  error  bound  over  the  measurements  points  i/m,m  =  1  Note  in  this  context  our 

hypothesis  reads 

Uk  =  {\Sk(um)\  <Bk,l<m<  M}.  (29) 

We  then  perform  the  statistical  test  of  Section  2  for  the  selected  confidence  level  7. 

If  we  accept  the  hypothesis  we  know  that  we  can  explain  the  difference  between  the  experimental  mea¬ 
surements  and  the  reduced  basis  prediction  in  terms  of  the  reduced  basis  approximation  error:  the  proposal 
k  should  thus  be  included  as  a  viable  (or  consistent)  candidate  for  k*  .  We  recall  that  we  can  not  conclude 
that  the  truth  model  is  indeed  valid  simply  because  we  do  not  reject  a  value  k:  we  do  not  control  the  Type 
II  error.  On  the  other  hand,  if  we  reject  the  hypothesis  then  we  know  with  confidence  greater  than  7  that 
the  proposed  truth  model  is  not  consistent  with  the  experimental  data.  In  particular,  the  discrepancy 

between  experiment  and  model  is  larger  than  can  be  attributed  to  the  error  incurred  in  the  reduced  basis 
approximation  of  the  truth: 

max  |sPS(r,m)  -  sh{vm,k) |  +  max  \sh{vm,  k)  -  sN(vm,k)\  >  max  |sPS(i/m)  -  sN(vm,k)\  >  maxAjv(i/m,K), 

I'm  I'm  Vm 

and  hence  since  max  \sh{i'm^  k)  —  sp]{vm,  k)|  <  maxAjv(^m,  k)  it  is  not  possible  that  sPS(^)  =  Sh{v,  k).  The 

I'm  Vm 

cause  for  the  discrepancy  must  thus  be  either  a  proposed  value  k  which  differs  from  the  actual  physical 
value  k*  or  some  physical  effects  that  are  simply  not  included  in  our  proposed  PDE  model.  Note  that  our 
hypothesis  can  not  distinguish  between  these  two  possibilities. 

We  may  now  consider  many  values  of  k.  (In  actual  practice,  we  would  presumably  only  search  in  some 
neighborhood  of  the  (or  a)  nonlinear  least-squares  “misfit”  minimizer  provided  by  (in  our  implementations) 
the  Levenberg-Marquardt  algorithm  [31].)  The  set 

VK  =  {kG  VK  |  HK  accepted  } 

constitutes  a  parameter  uncertainty  region,  or  more  optimistically  a  parameter  “possibility”  region:  VK 
includes  all  values  of  k  for  which  the  discrepancy  between  data  and  model  prediction  is  consistent  with 
the  reduced  basis  accuracy.  Note  that  VK  will  depend  on  the  accuracy  of  the  reduced  basis  approximation, 
the  values  of  7,  K,  M  and  I,  and  of  course  also  the  sensitivity  of  the  output  to  the  parameters.  The 
test  of  hypothesis  for  each  k  requires  evaluation  of  S]y(i'm,k),  1  <  m  <  M:  these  calculations  can  be 
readily  effected  thanks  to  the  extremely  low  cost  of  the  Online  reduced  basis  output  evaluation.  In  essence, 
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the  reduced  basis  output  approximation  permits  us  to  address  the  parameter  uncertainty  issue  —  direct 
calculation  of  the  truth  output  Sh{vm,  k),  1  <  m  <  M ,  for  many  k  in  VK  would  be  prohibitive  (certainly  in 
real-time);  the  reduced  basis  error  bound ,  in  conjunction  with  our  hypothesis  test,  permits  us  to  nevertheless 
provide  rigorous  conclusions  with  respect  to  the  truth  approximation. 

Of  notable  interest  is  the  case  in  which  the  parameter  possibility  region  is  empty:  for  all  k  £  VK  we  reject 
our  hypothesis.  (Note  we  can  no  longer  attribute  a  confidence  level  to  this  “joint”  outcome  since  our  Type  I 
error  analysis  considers  only  a  single  value  of  k.)  In  this  case  we  would  conclude  that  the  misfit  can  not  be 
attributed  to  reduced  basis  error  and  hence  the  physical  system  of  interest  can  not  be  adequately  represented 
by  the  proposed  truth  PDE  model  —  we  suppose  here  that  VK  is  sufficiently  large  (or  includes  all  plausible 
parameter  values)  such  that  we  can  safely  attribute  the  rejections  to  unmodeled  physics.  This  situation 
is  most  instructive  if  we  can  interpret  the  unmodeled  physics  as  changes  to  a  system  such  as  the  quasi¬ 
static  appearance  of  defects  or  geometric  variations  relative  to  an  initial  baseline,  or  perhaps  manufacturing 
departures  from  nominal  specifications;  if  our  hypothesis  is  rejected  for  all  k  then  we  can  conclude  that  the 
system  is  indeed  appreciably  different  from  our  reference  or  desired  configuation.  This  “detection”  capability 
can  be  very  effective  as  it  permits  us  to  exploit  a  relatively  low-dimensional  reduced  basis  parametrization 
x  VK  —  to  deduce  potentially  very  high-dimensional  (and  important)  deviations  —  for  example,  the 
appearance  of  defects  (e.g.  voids  or  cracks  or  tumors)  in  a  physical  system. 

Finally,  we  could  also  include  a  perturbative  term  on  the  right-hand  side  of  our  proposed  partial  differ¬ 
ential  equation  to  account  for  unmodeled  physics.  We  may  not  know  the  specific  details  of  the  perturbation 
(for  example,  spatial  fluctuations  in  the  thermal  conductivity  on  a  fine  scale),  but  based  on  plausible  as¬ 
sumptions  and  our  reduced  basis  stability  factors  we  may  be  able  to  bound  the  effect  of  the  perturbation 
on  the  output.  We  can  then  include  this  additional  “allowable  deviation”  in  our  choice  of  B  in  (21):  now 
B  will  be  the  sum  of  the  reduced  basis  error  and  unmodeled  physics  contributions.  We  do  not  consider  this 
extension  in  the  numerical  results  section  here,  but  it  is  a  possible  direction  for  more  robust  treatment  of 
(inevitable)  unmodeled  physics  in  future  work. 

4.  Numerical  Results 

We  illustrate  the  frequentistic  uncertainty  framework  developed  above  for  two  example  problems.  The 
first  example  is  a  transient  thermal  conduction  problem  with  uncertain  thermal  conductivity  parameters; 
the  second  example  is  a  Helmholtz  acoustics  problem.  A  libMesh-based  implementation  of  the  certified 
Reduced  Basis  method  [32,  33]  is  employed  to  generate  our  finite  element  and  RB  numerical  results. 
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Figure  1:  The  transient  thermal  conduction  problem.  The  origin  of  the  x  =  ( X\,X2 )  coordinate  system  is 
the  lower  left  corner. 

4-1.  Transient  Thermal  Conduction 
Proposed  PDE  Model 

We  consider  transient  thermal  conduction  in  the  domain  depicted  in  Figure  1.  We  consider  two  layers  of 
material  12 1  and  122  with  respective  parametrized  conductivities  K\  and  H2\  these  two  layers  are  affixed  to 
a  substrate  12o  with  specified  (normalized)  conductivity  kq  =  1.  We  impose  a  heat  flux  on  the  surface  and 
we  take  as  our  output  the  average  temperature  in  an  adjacent  “measurement  region”  as  a  function  of  time. 
We  now  present  the  the  mathematical  statement  of  the  proposed  parametrized  PDE  model.  Note  that  the 
proposed  parametrized  PDE  model  shall  include  both  the  equation  for  the  field  as  well  as  the  prescription 
for  the  output. 

We  consider  the  computational  domain  12  =  120Ul2iUl22,  where  r20  =  (2, 10)  x  (0, 10),  12i  =  (0, 1)  x  (0, 10), 
and  D2  =  (1,2)  x  (0,10).  We  define  Tin  =  {a:  €  <912  :  X\  =  0  and  X2  S  (4.5,  5.5)},  T00  =  {x  £  <912  :  X2  = 
0  or  X2  =  10  or  X\  =  10},  and  the  Sobolev  space  X  =  {v  £  U1(12)  :  v  =  0  on  T00}.  We  consider  the 
time  interval  t  £  [0,2/],  where  2/  =  5;  we  divide  the  time  interval  into  J  =  100  subintervals  of  equal 
length  A2  =  0.05  and  let  P  =  jAt,  0  <  j  <  J.  We  consider  the  two-dimensional  parameter  domain 
k  =  (ki,«2)  £  =  [0.2,  5]2.  Note  for  this  problem  v  from  Section  2  corresponds  to  time  (hence  =  0 

and  ^max  =  tf)  and  k  from  Section  3  corresponds,  in  fact,  to  n. 

We  first  consider  the  (semidiscrete  backward  Euler)  weak  form  of  the  governing  PDE:  For  n  = 
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(k i,  K2)  G  and  j  =  1, . . . ,  J,  find  k)  £  X  that  satisfies 


i(P,k)  —  u(P  1,«) 
At 


2 

+  /  Vu(tJ',  k)  •  Vv  +  /  KiVu(tJ ,  k)  •  Vv  =  /  w,  \/v  £  X.  (30) 

J  Oq  i 1  “  Fli  “  F in 


We  also  introduce  the  output  functional 


^0 


Vt;  G  X, 


(31) 


where  flme  =  (0,0.2)  x  (4.5,  5.5);  our  output  is  then  given  by  s(P ,  k)  =  £{u{P ,  k))  for  0  <  j  <  J. 

To  obtain  Xi^K  we  now  replace  X  in  (30)  with  a  truth  finite  element  space  Xh  C  X  with  dim(Afc,)  =  3009 
degrees  of  freedom:  For  k  =  («i,  n2)  £  VK  and  j  =  1, . . . ,  J,  Uh(P ,n)  G  Xh  satisfies 


uh(t3 ,  k)  -  uh{tJ  1,k) 
A  t 


+  /  Xuh(t3 ,k)  ■  Vu  + 


/  Oq 


2 

E 

i—l 


KiS7Uh(t3  ,  k)  •  Vc 


V«  G  Xh. 


(32) 


The  truth  output  is  then  given  by  Sh{t3 ,k)  =  £(uh(P ,k))  for  0  <  j  <  J. 

The  RB  approximation  un(P,k)  also  inherits  the  Galerkin  formulation  from  (30),  though  now  the  ap¬ 
proximation  space  is  the  RB  space  X^:  For  k  =  (k\,  H2)  £  T>K  and  j  =  1, . . . ,  J,  un(P  ,  k)  G  Xjy  satisfies 


un(P  ,  k)  —  un{P  1,k) 
A t 


1  +  f  Xun(P ,  k)  ■  Vv  + 
dn0  i=1 


KjXUN(P  ,  k)  •  Vc 


Vw  G  Aat.  (33) 


The  RB  output  is  given  by  Sn{P,k)  =  £(ui v(P,k))  for  0  <  j  <  J.  We  generate  Xn  via  the  Greedy 
algorithm  [21]:  we  require  dim(A,v)  =  40  in  order  to  satisfy  an  error  bound  tolerance  of  0.002  on  an  RB 
“training  set”  —  a  uniform  grid  of  900  training  parameters  in  T>K .  We  note  that,  for  given  n  G  T>K,  in  the 
Online  stage  —  relevant  to  many  queries  —  it  is  109  x  faster  to  evaluate  the  RB  ( N  =  40)  output  and 
output  error  bound,  k  — >  Sn(P ,  k),  An(P  ,  k),  than  to  directly  evaluate  the  truth,  k  —>  Sh(P ,  «)•  Recall  that 
\sh{P ,  k)  —  sjy(P ,  k)|  <  Ajy(P ,  k),  0  <  j  <  J,  k  £  An(P ,  k)  is  a  rigorous  error  bound. 

We  now  apply  our  statistcal  framework  to  this  proposed  parametrized  PDE  model.  We  consider  two  cases 
below.  In  the  first  case  we  assume  that  the  physical  system  (PS)  is  fully  modeled  by  the  proposed  PDE; 
we  perform  validation  hypothesis  tests  to  determine  which  values  of  the  conductivities  k  are  consistent  with 
(synthetic)  experimental  data.  In  the  second  case  we  consider  “unmodeled  physics” ;  we  introduce  a  physical 
system  (PS7)  which  does  not  correspond  to  the  proposed  parametrized  PDE  for  any  choice  of  parameters. 
Note  we  of  course  assume  that  we  are  unaware  of  the  “unmodeled  physics” :  the  validation  procedure  should 
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inform  us  that  the  physical  system  PS'  can  not  be  modeled  by  the  proposed  PDE. 


Fully  Modeled  Physics 

We  must  first  “construct”  our  Physical  System  PS.  Throughout  this  section  we  set  sPS(P)  =  Sh(V ,k*) 
for  the  reference  parameter  k*  =  (4, 1);  we  then  generate  synthetic  realizations  of  the  form  Vexp(t\  u>)  = 
sPS(P)  +  eG(P ,  w),0  <  j  <  J.  Here  the  Gaussian  process  G(t,uj)  is  synthesized  from  the  Karhunen- 
Loeve  expansion  (5)  with  correlation  length  either  A„  =  0  or  Au  >  0  (with  C(P,P  )  given  by  (6)  or  (7), 
respectively).  We  consider  the  proposed  PDE  (30)  so  that  the  RB  output  Sn  is  determined  by  (33),  and  our 
hypothesis  tests  are  based  upon  YK(P,  w)  =  [sPS(P)  —  Sn(P, «)]  +  eG(P,u).  We  consider  I  =  0  —  a  single 
Legendre  mode  —  unless  otherwise  indicated. 

We  first  demonstrate  the  effect  of  the  choice  of  N  on  the  possibility  region.  We  set  e  =  0.005  in  (25) 
and  generate  K  =  10  realizations  of  0  <  j  <  J,  for  a  Gaussian  process  G  which  is  temporally 

uncorrelated,  =  0;  we  consider  M  =  101  samples  in  time.  We  then  perform  hypothesis  tests  on  a  100  x  100 
uniform  grid  of  parameters  in  VK.  Figure  2  shows  the  (7  =  0.95)  possibility  region  for  N  =  10  and  N  =  40. 
The  reduction  in  B  which  results  from  an  increase  in  the  fidelity  of  the  RB  model  leads  to  a  much  sharper 
characterization  of  the  parametric  uncertainty  in  the  problem.  We  shall  henceforth  employ  only  the  N  =  40 
RB  space. 


(a)  (b) 

Figure  2:  Possibility  regions  for  e  =  0.005,  7  =  0.95,  K  =  10,  M  =  101,  1  =  0  and  (a)  N  =  10  and  (b) 
N  =  40;  k*  =  (4, 1)  is  indicated  in  each  plot  with  a  green  cross. 


We  next  consider  the  effect  of  K  and  M ,  still  with  uncorrelated  noise.  Figure  3a  is  a  repeat  of  Figure  2a 
to  provide  a  baseline  for  K  =  10  and  M  =  101.  In  Figure  3b  we  consider  K  =  2  and  M  =  101:  the 
reduction  in  I\  has  little  effect  here  because  the  noise  is  uncorrelated  —  the  large  M  compensates  for  the 
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small  K.  Figure  3c  presents  the  possibility  region  for  K  =  2  and  M  =  6,  and  in  this  case  the  possibility 
region  expands  signficantly  because  we  do  not  have  enough  data  samples.  In  short,  for  white  noise,  large 
M  leads  to  variance  reduction  because  our  calculation  of  the  “moment”  is  based  upon  M  independent 
samples. 
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Figure  3:  Possibility  regions  for  white  noise  with  e  =  0.005.  We  set  7  =  0.95,  /  =  0,  N  =  40  and  (a)  K  =  10, 
M  =  101;  (b)  K  =  2,  M  =  101;  (c)  K  =  2 ,M  =  6. 


We  repeat  the  above  study  of  the  effects  of  K  and  M  but  now  with  correlated  noise  defined  according  to 
(7)  with  =  1.  We  now  observe  in  this  correlated  noise  case  —  and  unlike  the  white  noise  case  —  that 
for  K  =  2,  Figure  4b,  the  possibility  region  is  significantly  larger  than  for  K  =  10,  Figure  4a.  Here  the  M 
temporal  samples  within  a  given  “experiment”  (or  realization)  are  correlated  and  therefore,  compared  to  the 
white  noise  case  of  Figure  3,  we  require  more  independent  realizations  (larger  K)  to  reduce  the  variance  of 
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Figure  4:  Possibility  regions  for  correlated  noise  (At  =  1)  with  e  =  0.005.  We  set  7  =  0.95,  I  =  0,  N  =  40 
and  (a)  K  =  10 ,M  =  101;  (b)  K  =  2,  M  =  101;  (c)  K  =  2 ,  M  =  6. 


Finally,  we  return  to  the  white  noise  case  and  demonstrate  the  effect  of  changes  in  the  choice  of  /.  We 
observe  from  Figures  5a  (/  =  0)  and  5b  (/  =  4)  that  for  I  =  4  we  obtain  a  much  sharper  characterization 
of  the  possibility  region  —  thanks  to  the  higher  fidelity  characterization  of  5(t).  However,  Figure  5c  shows 
that  for  1  =  8  we  again  increase  the  size  of  the  possibility  region:  the  culprit  is  the  denominator  of  (19); 
we  would  need  more  independent  experiments  (larger  K)  in  order  to  accurately  characterize  the  (/  =  8) 
9-variate  normally  distributed  data.  We  note  from  our  “best  result”  Figure  5b  that  there  is  greater  relative 
uncertainty  in  (less  sensitivity  of  the  output  to)  K2'  the  first  layer  enjoys  undivided  attention  for  the  shorter 
times  during  which  the  thermal  penetration  depth  has  not  yet  reached  the  second  layer.  Note  that  e  =  0.005 
corresponds  to  a  relative  output  noise  of  roughly  1%. 
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Figure  5:  Possibility  regions  for  white  noise  with  e  =  0.005.  We  set  7  =  0.95,  K  =  10,  M  =  101,  IV  =  40 
and  (a)  I  =  0;  (b)  I  =  4;  (c)  7  =  8. 


Unmodeled  Physics 

As  already  indicated,  in  our  unmodeled  physics  scenario  the  proposed  parametrized  PDE  is  described,  as 
before,  by  (30).  In  order  to  construct  our  unmodeled  physics  and  synthetic  “experimental”  data,  however, 
we  shall  require  a  new  PDE  model  which,  of  course,  must  include  the  unmodeled  physics  —  in  our  case  a 
crack:  this  new  PDE  model,  described  in  Figure  6,  shall  carry  a  superscript  “cracked”.  In  this  “cracked” 
PDE  model  we  introduce  a  crack  of  length  Lc  between  17 1  and  D2.  We  apply  zero  flux  boundary  conditions 
on  the  surface  of  the  crack;  the  weak  form  is  then  given  by  (30)  except  that  the  domain  17  is  replaced  by 
17'  (the  zero  flux  condition  is  natural).  Now  s“ack(t-7,  k*,  Lc)  shall  denote  the  truth  finite  element  output 

defined  in  the  same  way  as  (31)  except  now  the  integrals  are  performed  over  17'  —  for  the  “cracked” 
PDE  system  with  crack  length  Lc.  We  now  set  the  output  of  a  “cracked”  Physical  System,  PS',  to  be 
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Figure  6:  Transient  thermal  conduction  in  a  domain  which  contains  a  crack  of  length  Lc  between  Sli  and 

sPS  (P)  =  s^rack(t-J,  re*,  Lc);  we  again  set  k*  =  (4, 1).  The  experimental  realizations  are  now  generated  by 
Ve*p(P,uj)  =  sPS  (P)  +  eG(P ,  w),  0  <  j  <  J,  where  the  Gaussian  process  G(P,uj)  is  calculated  according 
to  (5).  We  again  consider  the  proposed  PDE  (30)  so  that  the  RB  output  sn  is  determined  by  (33),  hence 
YK(P,u})  =  [sPS  (P)  —  SAr(tJ,  k)]  +  eG(P,uj)  —  we  perform  model  validation  with  respect  to  the  baseline 
“idealized”  system  to  test  for  the  presence  of  pertubations  or  defects. 

Figure  7  shows  the  possibility  regions  with  cracks  of  length  Lc  =  1  and  Lc  =  1.5.  Figure  7a  illustrates 
that  the  Lc  =  1  crack  results  in  a  shift  in  the  possibility  region  so  that  it  no  longer  contains  «*;  we  note  that 
the  shift  is  predominantly  toward  lower  as  the  crack  results  in  a  reduction  in  heat  penetration  through 
fii.  Figure  7b  shows  that  with  Lc  =  1.5  we  obtain  an  empty  possibility  region;  hence  in  this  case,  with 
a  confidence  level  of  7  =  0.95,  we  have  deduced  the  presence  of  unmodeled  physics:  (32)  is  not  a  good 
model  for  PS'.  We  note  that  for  sufficiently  large  TV ,  I\ ,  and  I  it  would  also  be  possible  to  obtain  an  empty 
possibility  region  with  the  smaller  Lc  =  1  crack,  though  of  course  this  would  increase  the  computational  cost 
of  our  procedure. 
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Figure  7:  Possibility  regions  for  white  noise  with  e  =  0.005.  We  set  7  =  0.95,  K  =  10,  M  =  101,  N  =  40, 
I  =  4  and  (a)  Lc  =  1;  (b)  Lc  =  1.5. 


4-2.  Helmholtz  Acoustics 
Proposed  PDE  Model 

We  consider  a  Helmholtz  acoustics  impedance  model  in  the  domain  D  shown  in  Figure  8.  In  the  proposed 
parametrized  PDE  model,  [we  choose  as  parameters  the  frequency  of  the  (time-harmonic)  sound  field,  k,] 
and  the  complex  impedance  of  the  wall  Z  =  Zr  +  iZi  ( i  =  \/ — 1);  we  provide  a  sound  source  at  a  speaker 
diaphragm,  and  we  take  as  our  output  the  average  pressure  at  the  microphone.  Similar  to  the  previous 
example,  the  proposed  parametrized  PDE  model  shall  also  include  both  the  equation  for  the  field  and  the 
prescription  of  the  output. 

The  relevant  domain  boundaries  are  the  speaker  source  (Pin),  the  wall  of  unknown  impedance  (rw),  the 
truncated  exterior  (Prad),  and  finally  the  microphone  rmlc.  The  speaker  diaphragm  is  of  length  a  =  1.0,  the 
wall  is  of  length  L  =  12,  and  the  semicircular  truncated  boundary  is  of  radius  i?rad  =  6  (measured  relative 
to  the  midpoint  of  the  wall);  the  distance  from  the  diaphragm  to  the  wall  is  4,  and  rmlc  is  halfway  between 
the  diaphragm  and  the  wall.  Our  v  in  Section  2  is  now  k  with  "D^.  =  [1, 2],  and  n  of  Section  3  is  now  (Zr,  Z)) 
(impedance,  abbreviated  as  Z)  with  (. Zr,Zi )  £  Vz  =  [1,4]  x  [0,1]. 

We  first  consider  A/tl',K,  the  weak  form  of  the  governing  PDE:  for  ( k,Z )  £  T>k  x  2?z>  find  the  acoustic 
pressure  <j>(k,Z)  £  A'  that  satisfies 


V<j){k,Z)-Xv-  k2<f>(k,Z)v+ 


ik 


Zr  P  iZ; 


c t>(k ,  Z)v+ 


ik+ 


2 i?rad 


<t>{k,Z)l 


\/v  £  X, 
(34) 
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Figure  8:  Acoustics  impedance  problem. 


where  X  is  the  Hilbert  space  X  =  {v  G  H1( f2)}  (over  complex- valued  fields)  and  v  denotes  the  complex 
conjugate  of  v.  We  note  that  the  third  and  fourth  terms  of  the  left-hand  side,  and  the  right-hand  side  of  (34), 
represent  the  impedance  boundary  condition  on  the  wall  Tw,  the  first-order  radiation  boundary  condition  on 
the  truncated  exterior  boundary  rrad,  and  the  source  at  the  speaker  diaphragm  Fjn,  respectively.  We  then 
consider  the  output 


s(k,  Z)  =  £(<f>(k,  Z))  =  |j4tejRe  j  Jrniic  Z)  | ,  Vu  G  X,  (35) 

which  measures  the  real  part  (Re)  of  the  averaged  pressure  over  the  (microphone)  strip  rmlc. 

To  obtain  Xi'X'  we  now  replace  X  with  a  truth  finite  element  space  Xh  C  X  of  dimension  dim(A/j)  = 
14126  degrees  of  freedom:  For  Z  G  Vz  and  k  G  T>k,  <j>h(k,  Z)  satisfies 


X(j>h(k,Z)X7v -  k2(f>h(k,Z)v+ 


ik 


■iZ, 


<t>h(k,Z)v+ 


/prad 


ik+ 


2 i?rad 


<f>h(k,  Z)v 


Vu  G  Xh. 
(36) 


The  truth  output  is  then  given  by  Sh.(k,  Z)  =  £((f>h(k,  Z)). 

We  then  denote  the  RB  approximation  space  as  Xn  and  define  our  RB  approximation:  For  Z  G  T>z  and 
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k  £  T>k,  cj)N(k,Z)  satisfies 


J  S7(f>N(k,Z)  'Vv- J  k2<j>N(k,Z)v+J  ^  ^ ^  <f>h(k,  Z)v+ J  ^ik+^^j(j)N(k,  Z)v  =  -  j  ikv,  Wv  £  XN 

(37) 

our  RB  output  is  then  given  by  s^{k,Z)  =  £((j)N(k,  Z)).  We  generate  the  RB  space  XN  by  the  Greedy 
algorithm:  we  require  dim(A/v)  =  N  =  35  in  order  to  satisfy  an  error  bound  tolerance  of  10-5.  We  note 
that,  for  given  (fc,  Z)  £  T>k  x  Vz ,  in  the  Online  stage  -  relevant  to  many  queries  -  it  is  roughly  lOOOx  faster 
to  evaluate  the  RB  (N  =  35)  output  and  output  error  bound,  (fc,  Z)  — ►  SAr(fc,  Z),  A N(k,  Z),  than  to  directly 
evaluate  the  truth,  (k,  Z)  — >  Sh(k,Z).  Recall  that  \sh(k,Z)  —  S]y(k,Z)\  <  A N(k,Z),  ( k,Z )  £  Vk  x  Vz; 

A v(/c,  Z )  is  a  rigorous  error  bound. 

We  now  apply  our  statistical  framework  to  this  proposed  PDE  model.  We  also  consider  two  cases,  as  in 
the  previous  example:  a  “fully  modeled  physics”  case  in  which  we  assume  that  the  physical  system  (PS)  is 
fully  modeled  by  the  proposed  PDE  model  and  validation  hypothesis  tests  are  performed  to  determine  which 
value  of  the  impedance  is  consistent  with  (synthetic)  experiment  data;  and  a  “unmodeled  physics”  case  in 
which  a  physical  system  (PS^  is  introduced  which  does  not  correspond  to  the  proposed  parametrized  PDE 
for  any  choice  of  parameters.  Note  we  of  course  assume  that  we  are  unaware  of  the  “unmodeled  physics”:  the 
validation  procedure  should  inform  us  that  the  physical  system  (PS^  can  not  be  modeled  by  the  proposed 
PDE.  The  confidence  level  for  both  cases  is  7  =  0.95. 

In  what  follows,  we  shall  employ  a  similar  approach  as  in  the  previous  section  except  that  now  each 
realization  corresponds  to  measurements  at  M  different  frequencies  (rather  than  different  time  levels).  The 
interval  k  £  [1,2]  is  divided  into  J  =  7  subintervals  of  equal  length  Afc  =  1/7  and  we  set  k3  =  1  +  jAk, 

0  <  j  <  J;  measurements  are  performed  for  each  frequency  k3 ,  0  <  j  <  J. 

Fully  Modeled,  Physics 

We  first  “construct”  our  Physical  System  PS.  We  set  sPS(k3)  =  Sh{k3,Z*)  for  the  reference  parameter 
Z*  =  (. Z*,Z *);  we  then  generate  realizations  Vexp(k3,oo)  =  sFS(k3)  +  eG(k3,uj),  0  <  j  <  J,  synthesized 
from  the  Karhunen-Loeve  expansion  (5)  with  white  noise.  We  consider  the  proposed  PDE  (34)  so  that  the 
RB  output  Sn  is  determined  by  (37)  with  N  =  35,  and  our  hypothesis  tests  (for  any  given  candidate  Z)  are 
based  upon  Yz(k3 ,w)  =  \sPS(k3)  —  SN^k3-,  Z)]  +  eG(k3 ,oj).  For  all  the  results  below,  we  consider  7  =  0.95, 

K  =  8  realizations,  M  =  J  +  1  =  8  samples  in  frequency,  and  I  =  3  corresponding  to  /  +  1  =  4  Legendre 
modes.  Our  possibility  region  is  constructed  on  a  uniform  100  x  100  grid  of  parameters  in  T>z- 

In  Figure  9a  and  9b  we  present  the  possibility  regions  for  e  =  0.06  and  e  =  0.01,  respectively.  We  observe 
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that  the  reduction  in  noise  level  leads  to  a  much  sharper  characterization  of  the  parametric  uncertainty.  For 
this  “best”  choice  of  test  options  (K,  M,  I,  N )  we  obtain  a  rather  small  possibility  region  for  e  =  0.01.  Note 
that  e  =  0.01  corresponds  to  a  relative  output  noise  of  roughly  2%. 
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Figure  9:  Possibility  regions  for  (a)  e  =  0.06  and  (b)  e  =  0.01;  Z*  =  (3.5,  0.75)  is  indicated  in  each  plot  with 
a  green  cross. 


Unmodeled  Physics 

In  our  unmodeled  physics  scenario  the  proposed  parametrized  PDE  is  described,  as  before,  by  (34). 
In  order  to  construct  our  modeled  physics  and  synthetic  “experimental”  data,  we  shall  define  a  new  PDE 
which  includes  the  unmodeled  physics  -  in  our  case,  a  heterogeneous  wall:  this  new  PDE  model,  which 
corresponds  to  the  model  shown  in  Figure  10,  shall  carry  a  superscript  “2-part  wall.”  In  particular,  we 
now  assume  that  the  wall  consists  of  two  parts:  the  upper  part,  rWjt0p  with  real  infinite  impedance  (hard 
wall);  and  the  lower  part,  rw!bot,  with  impedance  (Zr,  Zi)  £  P>z-  The  “2-part  wall”  PDE  weak  form  is  thus 
the  same  as  in  (34),  except  now  the  integral  over  Tw  is  taken  only  over  the  lower  half  of  the  wall,  FWibot- 
We  now  set  sPS  (Ad)  =  part  wall(Ad,  Z*),  where  part  wall(Ad,  Z*)  denotes  the  “2-part  wall”  truth  finite 
element  output  for  Z*  =  (3.5,0.75).  Finally,  we  take  Fexp(Ad,  oj)  =  sPS  (Ad)  +  eG(Ad,  oj),  0  <  j  <  J, 
where  the  Gaussian  process  G(k,oj)  is  calculated  according  to  (5).  We  now  pursue  our  hypothesis  tests  with 

(Ad ,  oj)  =  [sPS  (Ad)  —  sjv(A:j';  Z)\  +  eG(A:J  ,  uj)  for  sn  evaluated  from  (37)  for  N  =  35:  we  perform  validation 
with  respect  to  the  baseline  homogeneous  wall  to  detect  the  presence  of  heterogeneities. 

Figure  11  shows  the  possibility  regions  with  noise  levels  e  =  0.06  and  e  =  0.01.  Figure  11  illustrates  that 
for  large  noise  the  difference  between  Fexp  and  Sn  results  in  a  shift  in  the  possibility  region  such  that  it 
no  longer  contains  Z*\  however,  for  small  noise,  we  obtain  an  empty  possibility  region.  In  the  former  case 
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the  test  (incorrectly)  reconciles  the  data  through  an  artificially  large  impedance.  In  the  latter  case  we  have 
deduced  the  presence  of  unmodeled  physics  -  a  wall  with  non-uniform  impedance. 


(a) 


(b) 


Figure  11:  Possibility  regions  for  (a)  e  =  0.06  and  (b)  e  =  0.01. 
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